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ABSTRACT 

In this paper, we implement a perturbative approach, first proposed by |Bouchet & Gispert| ( |I999l l, to estimate 
variation of spectral index of galactic polarized synchrotron emission, using linear combination of simulated 
Stokes Q polarization maps of selected frequency bands from WMAP and Planck observations on a region of 
sky dominated by the synchrotron Stokes Q signal. We bnd that, a brst order perturbative analysis recovers 
input spectral index map well. Along with the spectral index variation map our method provides a hxed 
reference index, /3oj, over the sky portion being analyzed. Using Monte Carlo simulations we bnd that, (/ios) = 
-2.84 ± 0.01, which matches very closely with position of a peak at Psip) = -2.85, of empirical probability 
density function of input synchrotron indices, obtained from the same sky region. Eor thermal dust, mean 
recovered spectral index, {fid) = 2.00 ±0.004, from simulations, matches very well with spatially hxed input 
thermal dust spectral index fid = 2.00. As accompanying results of the method we also reconstruct CMB, 
thermal dust and a synchrotron template component with hxed spectral indices over the entire sky region. 

We use full pixel-pixel noise covariance matrices of all frequency bands, estimated from the sky region being 
analyzed, to obtain reference spectral indices for synchrotron and thermal dust, spectral index variation map, 
CMB map, thermal dust and synchrotron template components. The perturbative technique as implemented 
in this work has the interesting property that it can build a model to describe the data with an arbitrary but 
enough degree of accuracy (and precession) as allowed by the data. We argue that, our method of reference 
spectral index determination, CMB map, thermal dust and synchrotron template component reconstruction is a 
maximum likelihood method. 

Subject headings: cosmic background radiation — cosmology: observations — diffuse radiation 


1. INTRODUCTION 

Observations of microwave signal by WMAP and Planck 
satellite missions, in the frequency range 23 GHz to 857 GHz, 
can be used to obtain a wealth of precise cosmological in- 
for mation (e.g, precise values of basic cosmological parame¬ 
ters (|Hinshaw et al.|2013[|Bennett efak 201 3[ [Planck Collab^ 
oration et al.|2015c|), topology of space (Planck Collaboration 
et al.|2015t] ), laws and nature of physical conditions existing 
at different epoch marke d by different phases of dynamical 
evolu t ion of s pace-time ( jPlanck Collaboration et al.|2015g|e| 
|2014| |2015d| l) as well as precise information about various 
astrophysical mechanisms operating since re cent times that 
cause diffuse emissions inside Milky Way ([Bennett et al.j 
[2003 1 [Planck Collaboration et al.|2015h[[2013| [Dobler|2012| l. 

However, extracting these information accurately with preci¬ 
sion requires a method that separates different components 
from their observed mixture using a model that is flexible 
enough to describe the data with an arbitrary but desired de¬ 
gree of minute closeness as permissible by data. A com¬ 
pletely accurate model may produce ambiguous results if the 
data is simply insufficient to constrain all model parame¬ 
ters correctly e.g., consider the astrophysical problem of re¬ 
construction various galactic emission components which is 
a greatly complex task due to inherent degeneracy existing 
in the morphological pattern of these emissions, degeneracy 
arising due to similarities in physical models of multiple emis¬ 
sions or due to other reasons. Such uncertainties in recon¬ 
structed foreground components may cause cosmological in- 
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ferences drawn from the analysis of reconstructe d Cosmic 
Microwave Background (CMB) maps to be biased ( [Tegmark] 
[et al.|2000|l if these uncertainties are not known or if they are 
not properly taken into account. 

One of the prime source of uncertainties against prop¬ 
erly extracting the primordial signal and other galactic fore¬ 
ground components from observed maps is intrinsic varia- 
tion of spectral index of synchrotron emission ( |Tegmark|1998| l 
over the sky due to variations of number density of cosmic 
ray electrons in our galaxy. In this paper we implement a 
novel method on simulated WMAP and Planck Stokes Q fre¬ 
quency maps to perturbatively model the synchrotron spectral 
in dex variation over the sky f ollowing the approach proposed 
by [Bouchet & Gispert[ ( [1999l ). The perturbative nature of the 
method is very desirable since it allows one to build a model 
with just enough accuracy that is still meaningful within the 
available data. As associated results of the method we recon¬ 
struct best-fit map for CMB and templates for thermal dust 
and synchrotron emission. 

Usefulness of the method described in this paper can also be 
seen from a different perspective. There exist several methods 
in the literature for CMB or (and) foreground component re¬ 
construction which can be grouped broadly in following three 
categories: 


1. The so-called blind, internal-linear-combination ap- 



proach, ((Tegmark & Efstathiou 

1996 Maino et al. 


2002 

Tegmark et al.[ [20031 

Delabrouille et al. 


2009 

[Saha[[2011||Bennett[[2012[[ 

^lanck Collaboration 

1 et al. 

|2015a]l), where one removes foregrounds based 


upon the fact that CMB follows blackbody spectrum 
whereas, other components do not. In principle, this 
method can also be used to produc e a cleaned map of 
any of the foreground components ([Remazeilles et all] 














































2 


(|201 1|)), provided its spectral index is a priori known 
and is constant over the sky portion being analyzed. 
The method by dehnition is a non-Maximal Likelihood 
(ML) approach since it incorporates minimum informa¬ 
tion about the physical components and no information 
about detector noise. 

2. The second class ([Eriksen et al.| ( f2008a|b| l; [Planck Col-| 
laboration et aLl ( |2015b|a| l) is that one where one uses 
as much prior information as possible about all compo¬ 
nents, along with the detector noise model, to perform 
a ML analysis of the data. The outcome of this method 
are ML estimates of all components. 

3. The third one is an intermediate approach of the first 
and second categories. Although, here one incor¬ 
porates informations about various components avail¬ 
able from prior observations, (e.g.. Maximum Entropy 


Method, dHobson et al.l 19981 IBennett et al. 

2003HHin-| 

shaw et al. 2007 |Gold et al. 2009 

2011 

^ and the 

Wiener filtering approach, (|Bunn et a 

1.1 1994 

Tegmark 

& Efstathiou|1996HBouchet et al.|1999|l), it is a non ML 


approach. 


All the above methods have merit and demerit of their own. 
ML approach has advantages that it is the optimal method in 
the sense that it retains all information of the data and provides 
maximum precision on the derived parameters. ML method 
may, however, be biased if the model to be used for the data 
is not accurately known, or, if there are degeneracies between 
the parameters. The non-ML approaches are advantageous 
since data is processed with the minimal assumptions. How¬ 
ever, it is a sub-optimal method and can be biased. One, there¬ 
fore, desires to ask a question can we develop a method which 
requires minimal assumptions in terms of data model (e.g., no 
requirement of using templates) and, at the same time, which 
is a ML approach with respect to at the least almost all of 
its parameters? So far, there has not been any attempt in 
the literature to answer this question. The novel perturbative 
approach for component separation and for determination of 
variation of spectral index of synchrotron emission followed 
in this paper, at one-side, is a linear-combination method, on 
the other side, is a ML approach with respect to all of its com¬ 
ponent templates and best-ht spectral indices for synchrotron 
and thermal dust. The method incorporates information about 
spectral properties of components and detector noise model, 
thus preserving properties of first and second categories de¬ 
scribed above. 

There have many attempts in the lite rature to measure syn- 
chrotron spectral index of our Galaxy. jReich & Reich ( 


find map of synchrotron spectral index between 420 M. 

1420 MHz using radio continuum su rveys. Following the 
method of TT plot|Davies et al.|(| 1996)1 find synchrotron tem¬ 
perature spectral index lies in tne range -2.8 to -3.2 from a 
strip a t ^ = 40° latitude sect ions of 408 MHz and 1420 MHz 
maps. [PTatania et ah ( |1998 1 find synchrotron s pectral index 
betwee n 420 MHz and 7.5 GHz as -2.76 ±0.11. [Bennett etal.l 
( |2003| ) find a flatter synchrotron spectral index, ~ -2.5, near 
the galactic plane and a steeper index, ^ -2.9 towards the 


higher galactic latitude. Using WMAP nine-year data Fuske- 
land et al. (2014|l estimate polarized synchrotron spectral in- 


^ See also all but first one of these papers for analysis following methods 
belonging to first and second categories. 


dex towards galactic plane as -2.98 ±0.01, and a steeper spec¬ 
tral index, -3.12 ±0.04, towards the higher galactic latitude. 

We present our paper as follows. We describe the prob¬ 
lem in Section]^ In Sectionwe describe the perturbative 
approach used to model the data. We discuss basic formal¬ 
ism of this work in Section]^ Here we discuss 1) the power 
minimization algorithm (based upon the model described in 
Section and how the power minimization algorithm leads 
to an equation for weights that we use for linear combination 
of input Stokes Q maps, later in Section 2) the difference 
map at each frequency bands, and the noise covariance matrix 
thereof, 3) the data-model mismatch statistics and 4) rela¬ 
tionship of our method with ML approaches. In Sectionl^we 
describe the method used in this paper. We present results in 
Section]^ Finally, we conclude in Section]^ 

2. BASIC PROBLEM 

Let us assume we have observations of CMB and fore¬ 
ground polarizations at n* different frequency bands. Each 
frequency band contains contribution from n different physi¬ 
cal components, namely, CMB and emissions from « -1 fore¬ 
ground components - each with a given spectral index all over 
the sky - plus detector noise If foreground component j 
has a hxed spectral index, jdj, all over the sky, net polarized 
Stokes Q signal, Qfp), at a pixel p, due to CMB, Qdp), all 
foregrounds, Qqj, and detector noise, Q'l(p) at a frequency u/ 
in thermodynamic temperature unit is given by, 

\ 

Qi(p) = ec(F)+Vfl/ — Qoj(p)+Q'l(p). (1) 

where a, = a(z/,) denotes conversion factor from antenna to 
thermodynamic temperature unit at frequency z/, for all fore¬ 
ground components. CMB anisotropy at any given pixel p 
is independent on frequency iz, due to its blackbody nature. 
Each foreground template Qoj is defined with reference to fre¬ 
quency 1 ^ 0 j, which may be different for different foreground 
components. We note that, in general, observed maps for dif¬ 
ferent frequencies are convolved by different instrument re¬ 
sponse functions. However, we implicitly assume that fre¬ 
quency maps, Qiip), in Eqnfl] already have been brought to a 
common resolution by propmy deconvolving first by respec¬ 
tive instrument response functions and then convolving again 
by the common response function. Since all frequency maps 
of our work are convolved by a common Gaussian beam func¬ 
tion of EWHM = 20° we omit reference to beam function in 
all equations of this paper containing frequency maps, CMB 
or foreground templates. 

Eor polarized signal there are only two foreground com¬ 
ponents at microwave frequencies, namely, synchrotron and 
thermal dust. Moreover, the spectral index for synchrotron 
emission varies with sky position due to variation of relativis¬ 
tic electron density over the sky. In the presence of variation 
of spectral index for synchrotron, but fixed spectral index for 
thermal dust component Eqn[T]reduces to, 

/ s A(/» 

Qiip) = Qcip) + ai — Qosip) 

W) ^2) 

+ ai(—^ Qodip) + Q’!(P)^ 

\VQd J 

^ Detector noise does not count as a physical component. 
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where /3i(p) represents the synchrotron spectral index at a 
pixel p and f3d is the global thermal dust spectral index. The 
motivation for our work is to obtain a reliable estimate of 
Psip) using a perturbative analysis, without assuming any 
prior template models for foregrounds or CMB. 


3. PERTURBATIVE DATA MODEL 

Let US consider a given region of sky. Synchrotron spec¬ 
tral index at a pixel p can be written as, Psip) = /3o.s + ^/3.s(p)^ 
where /3 oj represents some reference spectral index for the 
region and A^sip) represents actual variation of spectral in¬ 
dex over /3oi inside the region. Using the direction dependent 
spectral in dex, (3s(p), in the second term of Eqn|^ one can 
show that ( |Bouchet & Gispert|1999| l, 


-) 


0s(p) 





^l-i-A/3^(p)ln 



(3) 


where the second line follows from Taylor expansion of 
the left hand side of above equation up to first order in 
Ajdslp)/Pos, assuming A/3^(p)//3os < 1. Using Eqn. in 
Eqn. we obtain. 


Qi(p) = Qc(p) + ai 



A, 

QQsip) + ai 



X 2os(p)A/3j(p)-i-a,- 



QM + aiip) 



Apart from CMB and detector noise, our data model given 
by above equation, can now be interpreted nicely to con¬ 
sist of a total of three foreground emission components, 
namely, synchrotron with a constant spectral index /3o.5, a sec¬ 
ond component 2oj(p)A/3j(p), with its frequency variation 

a, (^) In^^^ and hnally thermal dust component. As 

we can see from above equation our model renders the varia¬ 
tion of spectral index of synchrotron component in terms of a 
new foreground component with rigid frequency scaling. We 
define a set of templates, Tj, where j G {1,2,3,4} following. 


Ti{p) = Qc(p) 

Tiip) = Qos(p) 

Tiip) = Qos{p)Al3,(p) 
Tdp) = Qdip) ■ 


shape vector are given by|3 

s/= (^) a(vi), (7) 

where fdj represents spatially fixed spectral index of /* (fore¬ 
ground) component. Eollowing method as described in Sec- 
tion[^we reconstruct each template dehned by Eqn[^ Due to 
presence of detector noise, the reconstructed templates, Tj(p), 
are, however, different from Tj{p). Once all template compo¬ 
nents are reconstructed, synchrotron spectral index variation 
map is obtained following, 

A^{p) = %{p)/Up). (8) 

4. FORMALISM 

4.1. Component Reconstruction and Weights 

From a mixture of n physical components present in dif¬ 
ferent frequency bands, our aim is to isolate best-fit map of 
each one of them by removing rest. This can be achieved 
by reconstructing any one component hrst and then repeating 
the method for all other components. Let us denote the com¬ 
ponent to be reconstructed by a Greek letter subscript, e.g., a 
and other components to be removed by any Roman letter sub¬ 
script taking values from the difference set{l,2,...,n}-{Q!}. 
We denote shape vector for the component to be reconstructed 
by [s]q while each of the rest n -1 components has shape vec¬ 
tor [s];. To reconstruct the cleaned template, Ta{p), for com¬ 
ponent a, we form a linear combination of available maps 
following, 

lib 

fc.(p) = J2<Q‘(P^^ (9) 

!=1 

where the linear weights, w‘^, for a component is obtained by 
minimizing the total pixel-pixel variance, cr^, of the cleaned 
template, Taip). Pixel variance of of this template is defined 
as, 

(>»> 

where Afpix denotes total number of pixels available for analy¬ 
sis in a single frequency map after any suitable mask has been 
applied . Using matrix notation, above equation can easily be 
written as, 

^^ = [W]JV][W]^, (11) 


With these definitions, total polarized emission of Eqn.[^ at a 
frequency i/, is then given by, 

4 

Qfp) = J2Tj(p)si + Q’;{p), (6) 

,/=i 

where we have denoted emission from component at chan¬ 
nel i (in thermodynamic temperature unit) by Tj(p)sj. Tfp) 
denotes template for component based on reference fre¬ 
quency VQj and sf for a given j, denotes elements of so-called 
shape vector, [s]/, for component j. For foregrounds with 
emissions following usual power law behavior elements of the 


where, [W]a is a 1 x weight matrix, (wj^,w^,....,wj!,'’), 
for component a. [V] represents, nj, x nt covariance matrix, 
with elements V,,' representing covariance betwee n frequency 
maps Q i(p) and Qi'ip'). As described in details in Saha et ak] 
(20081 the choice of weights which minimizes cr^ is given by. 


[W]„ = 


[s]a[V]^ 
[s]a[V]ns]?; ’ 


( 12 ) 


where [V]'^ repre s ents Moore-P enrose generalized in¬ 
verse (|M oore||192Q| |Penrose||1955| > of [V]. Our aim is to 


^ For CMB elements of shape vector are unity for all frequency bands and 
choice of uqc is irrelevant. 
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obtain an analytical expression for weights that minimizes 
variances due to physical components in the limiting situa¬ 
tion when detector noise is negligible]^ Using variance and 
covariance of different components we can write, 

[V] = ai[s]^[s]„ + [r]^[s]„ + [s]^[r]„ + [F], (13) 

where [r]a is a 1 x n* matrix describing random chance cor¬ 
relation between the component to be reconstructed and all of 
the other components at each of the tii, frequency bands. [F] 
describes the rih x rii, covariance matrix of n - 1 components, 
hence it has a rank n -1. Above equation can be written as, 

[V] = a2[s]^[s]„ + [A]i, (14) 


Hence the total noise, Q"(p), due to all recovered components 
at frequency v,, is given by, 

n n rib 

Q'lip) = H 

k=l k=\ j=l 

(18) 

J=1 

where we have defined, 

n 

w\ = ^s’lw{. (19) 

k=\ 


where [A]i = [r]^[s]„ + [A ]2 and [A ]2 = [s]^[r]„ + [F]. Using 
above eq uation we can apply generalized Sherm ann-Morisson 
formula (|Baksalary et al.|(|2003|l;|Meyer|(|l973|)) successively 
to obtain both numerator and denominator of Eqn[T2](e.g., as 
in Saha et al. ( 2008| l). After a long algebra and negle^cting all 
terms containing [rjo,, Eqn[T^reduce to 


We form a difference map at each frequency, Vi, by subtracting 
the recovered and input frequency map at the same frequency, 
Vi. Using Eqn[^noise in the difference map is given by, 

rib 

Q"{p)-mp )=(1 -wm(p)-J2^iQ^ip) ■ m 

j¥i 


[W], 


[s]a([I]-[F][Fr) 

[s]a([I]-[F][F]+)[s]? 


(15) 


The term [F][F]'^ can be interpreted as the projector on the 
column space of [F]. Hence ([I] - [F][F]'^) is the orthogo¬ 
nal projector on the null space of [F]. We note that, [s], G 
C([F]),fori = {l,2,3,Hence, [F][Fr[s]f = [s]f 
i.e., [W]Q[s]f = 0, however, since [s]o ^ C([F]) we have 
[W]q[s]^ = 1, which is the condition for reconstruction of the 
component under consideration. Eqn. [TS] is not directly us¬ 
able by us since a priori we do not have any knowledge about 
the templates to be reconstructed from data and hence about 
the covariance matrix [F]. We, however, see that the crucial 
properties of [W]q,, viz., [W]a[s]f = 0 and [W]q[s]^ = 1 that 
must be satisfied for reconstruction of a component remain 
unchanged if we replace [F] in Eqnby a new rib x n-1 
matrix, [G], whose column vectors are given by each of the 
shape vectors of the n-1 components which are to be removed 
from the data. We therefore compute weights replacing [F] by 
[G] in Eqn[^ 


4.2. Difference Map and Noise Covariance 

Once all components have been reconstructed, using 
weights following Eqn for the complete set of possible 
shape vectors of all components, we can form reconstructed 
frequency map, Qiip), at each frequency, v,, following, 

n 

Qi{p) = Y,fj(p)s\. (16) 

7=1 

Due to underlying linearity in our method, one can clearly see 
that noise, Ql{p,Vi), in the reconstructed template at fre¬ 
quency Vi, has contribution from noise of all input frequency 
maps. Specifically, 


lib 

Ql{p,Vi) = s1Y,w{Q]{p). (17) 

7=1 

® The assumption of limiting detector noise is justified for our current 
work, since we use low resolution maps for our analysis. 


Assuming noise of different frequency bands have zero mean 
and noise properties of different detectors are uncorrelated 
with one another, noise covariance of the difference map 
is given by. 


Nf,! = (1 - mp)Qlip '))+E ('^/) {Q]ip)Q!]ip')) 

j^i 

nb , (21) 

= (i-w')2a;,,+E('^0 Kp'^ 


j¥i 


where is the noise covariance matrix of the i‘^{f^) input 
frequency map. 


4.3. Data-Model Mismatch Statistic 

Using the noise covariance matrix obtained in the previous 
section we define mismatch between the data and model by 
means of usual statistic following, 

^ E {Qi(P)-Qi(P)) {Qiip')-Qi(p')) (22) 

UP.P' 


where represents Moore-Penrose generalized inverse 


of Npp, and N = nbX A/j,ix denotes total number of surviving 
pixels counting all frequency bands after any suitable mask 
has been applied_on the maps. Since, Qi{p) depends upon 
Pd through Eqn 


16 


and depends upon /3oi, /3d and fre¬ 
quency band noise covariance matrices through Eqn we 
have = X^(/3o.'i,l3d)- Clearly, /3 oj and /3d enter as the min¬ 
imizing variables for the x^ statistic. Our best fit values for 
these spectral indices are the ones which together minimizes 
the x^ given by above equation. 


4.4. Relationship with ML approaches 


Eqn 22 shows that the best ht spectral indices are ML es¬ 
timates for a known detector noise model provided we have 
also chosen uniform prior on the indices. The ML nature of 
these indices directly implies that all the reconstructed tem¬ 
plates, Tj{p), are also ML estimates. We, however, note that 
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TABLE 1 
Frequency maps 


Frequency, u 
(GHz) 

WMAP, Planck 

Name 

a(u) 

aoipK) 

Ref" 

23 

WMAP 

K 

1.014 

1435.0 

30 

Planck, LFI 

A 

1.023 

36.56 

33 

WMAP 

Ka 

1.029 

1472.0 

41 

WMAP 

Q(l,2) 

1.044 

2254.0, 2140.0 

44 

Planck, LFI 

B 

1.051 

37.03 

61 

WMAP 

V(l,2) 

1.099 

3324.0, 2958.0 

71 

Planck, LFI 

C 

1.136 

37.11 

100 

Planck, HFI 

D 

1.283 

15.83 

143 

Planck, HFI 

E 

1.646 

11.80 

217 

Planck, HFI 

F 

2.961 

19.39 


Note. — List of IVMAP and Planck frequency maps used 
in the analysis. Wherever multiple differencing assemblies are 
available for WMAP polarized noise level for each differencing 
assemblies is also mentioned in orde r. Planck noise lev els are 
consistent with the Planck blue book ( |consortium|2005^ values 
and scaled to Wide = 512. For WMAP actual noise level per pixel 
at Wide =512 is given by a{p) = ao/^Nobs(p), where Wbs(p) 
denotes effective number of observations at a pixel p for given 
DA map. For Planck bands noise level per pixel at Wide = 512 is 
assumed to be uniform with (j(p) = (tq as the specified values at 
Wide = 512. 

“ Stokes (Q, U) polarization noise level in pK thermodynamic 
temperature unit at Wide =512. 

our reconstructed spectral index map is not a ML estimate 
since we construct it following Eqn|^ by taking ratio of two 
templates, without directly estimating it from the fit. 

5. METHODOLOGY 



0.02 


0.015 


- 0.01 

- 0.005 

- 0 


Fig. 1.— Full pixel-pixel noise covariance matrix obtained from A^nside = 8, 
20° smoothed maps over surviving pixels of pmlO mask for V band. 

noisy reconstruction of spectral index map, we remove those 
pixels from all 10 input frequency maps which contain weak 
synchrotron Stokes Q polarization signal. For this purpose, 
based upon the smoothed Stokes Q synchrotron template at K 
band, we make a mask which excise all pixels with absolute 
values less than 10pK in antenna temperature unit. We call 
this mask pmlO mask and apply this mask to all frequency 
maps for reconstruction of spectral index map[^ 


5.1. Input Data 

5.1.1. Frequency bands 

We include simulated Stokes Q polarization maps of 
WMAP K, Ka, Q, V bands, all three Planck LFl bands, and 
three lowest Planck HFI bands. We do not include WMAP W 
band in our analysis since it contains somewhat larger noise 
level compared to other WMAP bands. We do not include two 
highest frequency Planck HFI maps since they are strongly 
dominated by thermal dust component. Moreover, with in¬ 
clusion of HFI 353 GHz map in our analysis, values be¬ 
come increasingly more sensitive to thermal dust model, giv¬ 
ing lesser weightage to synchrotron component. Hence we 
exclude this frequency band also from our analysis. We how¬ 
ever verify that conclusions drawn in this paper from spectral 
index analysis remain unchanged with inclusion of HFI 353 
GHz map. A detailed specifications of the frequency bands 
used in this work is given in Table [1] 

Since pixel values of spectral index variation map, A/3s(p), 
are small compared to the foreground emission, any signifi¬ 
cant residual noise in the reconstructed template maps results 
in its potentially noisy reconstruction following Fqn|^ Ow¬ 
ing to delicate noise sensitivity of reconstructed spectral index 
map, we chose to minimize detector noise at the beginning of 
analysis, by first downgrading all Stokes Q frequency maps 
at Vside = 8 and then smoothing each one by a Gaussian beam 
function of FWHM = 20°. 

5.1.2. Mask 

We see from Fqn^ Aj3s(p) map may contain noisy pixels 
if corresponding pixds in the reconstructed synchrotron tem¬ 
plate, Ta{p), (or T-iip)), are also noise dominated. To avoid 


5.1.3. Physical Components 

We use WMAP ‘base fit’ MCMC (Markhov Chain Monte 
Carlo) synchrotron and thermal dust polarization maps for 
our simulations. These maps are provided at Aside = 64 at 
LAMBDA website. As mentioned at LAMBDA website, for 
some of the pixels of these maps, MCMC fitting resulted in 
errors. For both Q and U polarized synchrotron map at K 
band, using WMAP's bit-coded error flag file, we first identify 
all pixels for which WMAP MCMC fit resulted in error. The 
process identifies a total of 1315 pixels at Aside = 64. Pixel 
value for each of these bad pixels is then replaced by the aver¬ 
age pixel value of its nearest neighbors. If a bad pixel has one 
or more bad neighbours we average only over the good ones. 
Following similar procedure we obtain polarized thermal dust 
template at W band at Aside = 64. We downgrade each of these 
two templates to Aside = 8. We generate polarized CMB maps 
compatible to WMAP LCDM spectrum directly at Aside = 8. 
We smooth both foreground templates and CMB maps by the 
Gaussian beam of FWHM = 20°. 

5.1.4. Spectral Index Map 

For spectral index we use WMAP ‘best-fit’ MCMC spec¬ 
tral index map. This map is provided by the WMAP science 
team at Aside = 64. We downgrade this map at Aside = 8 and 
subsequently smooth by a Gaussian beam of FWHM = 20°. 
This smoothed spectral index map defines our primary spec¬ 
tral index map, /3j(p), for this work with which we compare 
the reconstructed spectral index map obtained by our method. 

^ We apply pmlO mask only for reconstruction of spectral index map. We 
present resu lts for reconstruction of CMB and foreground components (e.g., 
Section[63) over entire sky. 
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5.1.5. Detector Noise Maps and Covariance Matrix 

Generating noise maps and thus modelling the noise prop¬ 
erties in the low resolution maps are somewhat involved than 
simple smoothing operations that were needed to be per¬ 
formed on for eground templates and CMB maps at A^side = 8 in 
Section l5.1.3l We take into account WMAP's scan modulated 
detector noise pattern by using Nohsip) information associated 
with each nine-year Differencing Assembly (DA) map avail¬ 
able at LAMBDA website at Wide = 512. We first compute 
noise variance at each pixel following a^ip) = ctQ/Nohsip) at 
Wide = 512 for all DA maps between K to V band. Since for 
each of Q and V bands two DA maps are available we sim¬ 
ply average the variances of the corresponding DA’s at each 
pixel, to obtain noise variance map corresponding to average 
DA map within the frequency bands. For Planck LFI and 
HFI maps we use simple uniform white noise model. Pixel 
white noise levels, a{p), for Planck maps are obtained from 
Planck blue book specifications, abb, which represents noise 
in the beam area of the corresponding map in unit of fiK/K. 
Pixel noise levels for Planck maps at Wide then given by, 
a{p) = abbTcMBBfwhm/ A0, where Bf„hm represents FWHM of 
Planck beams in arcmin, A0 = 6.9', represents pixel size of 
Wide = 512 map and Tcmb = 2J3K. 

In principle, one can use all the noise variance maps at 
Wide = 512 corresponding to all frequency bands of Table 
to obtain representations of noise maps at Wside = 512 and 
then downgrade them to Wide = 8 for use in our simula¬ 
tion. However, one needs to simulate many such maps at 
Wide = 512, since we need to obtain pixel pixel noise covari¬ 
ance at Wide = 8 by downgrading the high resolution noise 
maps. Since generating noise maps at Wide = 512 and sub¬ 
sequently downgrading them many times are computationally 
very expensive we device a completely equivalent approach 
which is computationally significantly fast compared to the 
former method. Instead of downgrading Wide = 512 noise 
maps we obtain noise variance maps at Wide = 8 by simply 
averaging over noise variance of all Wide = 512 pixels lying 
inside each Wide = 8 pixel. Such a method is always feasible 
for our method and does not introduce any pixel-pixel corre¬ 
lation at Wide = 8 since noise in the original Wide = 512 maps 
are uncorrelated from pixel to pixel. Using these noise vari¬ 
ance maps we simulate noise maps at Wide = 8. These noise 
maps are then smoothed by Gaussian beam of FWHM = 20°, 
to generate realistic noise maps at Wide = 8 for each frequency 
band. We obtain 10000 Monte-Carlo simulations of smoothed 
noise maps at Wide = 8. Since the smoothing operations intro¬ 
duces pixel to pixel noise covariance we obtain full pixel to 
pixel noise covariance matrix over surviving pixels (after ap¬ 
plication of pmlO mask) at Wide = 8 corresponding to each 
frequency band. We have shown noise covariance matrix for 
V band masked map, in Figure Non-diagonal nature of this 
matrix indicates non-trivial noise properties in the maps. 

Using the Wide = 8 smoothed CMB Stokes Q polarized map, 
Qcip), synchrotron and thermal dust templates (Qos and Qod 
respectively), spectral index map, /3s(p) and detector noise 
map, Q”(p) (i = l,2,...,n*) as obtained above we generate 
simulated frequency maps for WMAP and Planck frequency 
bands, using Eqn S and knowing the values of a,. We use 
j3d = 2.0, vqs = 23 GHz and r^od = 94 GHz in this work. 

5.2. Analysis 

Since spectral indices of neither synchrotron, nor thermal 
dust components are known a priori, we perform reconstruc- 



Fig. 2.— Top panel. Mean synchrotron reconstructed spectral index map 
obtained from 1000 simulations of this work. Middle panel. The standai'd 
deviation of spectral index map. Bottom Panel Difference of top panel and 
input spectral index map. 

tion of all components defined by Eqn for a suitable range 
of both of these indices. We vary Pqs between -2.65 to -3.05 
with a step of -0.01. Eor each of these values is var¬ 
ied between 1.80 and 2.20 with a step size 0.01. Eor each 
pair {f5os,/3d), such generated, we compute shape vectors for 
all foreground template components. Combining these with 
the shape vector of CMB we estimate a set of weights us¬ 
ing Eqn [T^ to reconstruct all template components includ¬ 
ing CMBM^or each pair (/3os,Pd), we form reconstructed fre¬ 
quency maps following Eqn using all the reconstructed 
template components. We then form the difference map at 
each frequency by subtracting reconstructed frequency map 
from the corresponding input frequency map. The noise co- 
variance of the i'^ difference map over the surviving pixels 
for a given pair of indices is given by the matrix, Npp/, de¬ 
fined by Eqn|^ Once the difference map is formed we ob¬ 
tain x^(fios,Pd) for all indices following Eqn|^ Our best fit 
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values for indices are the ones which minimizes over all 
indices. Using these values we obtain spectral index variation 
map following Eqn[^ We validate our method by performing 
1000 Monte Carlo simulations of entire procedure described 
above. We note that since CMB follows blackbody distribu¬ 
tion its shape vector remain identical for all foreground in¬ 
dices and for all simulations. Since the weights as given by 
Eqn depend only upon foreground and CMB shape vec¬ 
tors (wthout any dependence on random detector noise which 
vary from simulation to simulation) for reconstruction pur¬ 
pose of all components, we calculate these weights only once 
for the entire range of spectral indices and store them on the 
disk. We then use these weights to reconstruct all components 
for all Monte-Carlo simulations for all indices. Moreover, we 
also precompute Moore-Penrose generalized inverses of de¬ 
tector noise covariance matrices of all difference maps for 
all pairs of spectral indices (/3os,(3d) and store them on the 
disk. In effect a Singular Value Decomposition algorithm with 
cutoff 1.0 X 10“^ was used to implement Moore-Penrose in¬ 
verses. 


6. RESULTS 

6.1. Spectral Index Map 

Eor each simulation described in Section we construct 
spectral index variation map, A/3s(p), over the reference in¬ 
dex, /3 oj, following Eqn|^ The results for reconstruction of 
spectral index map are shown in Eigure]^ The colored area of 
Mollweide projection of all the three panels shows the region 
on the sky that survives after applying pmlO mask. The top 
panel shows mean reconstructed synchrotron spectral index 

map, defined as, = (/3oi,min) + 

where /3o.!,min denotes value of /3os (e.g., see Eqn corre¬ 
sponding to Xmin the average is performed over all 1000 
simulations. The middle panel is simple standard deviation 
map for any one of the reconstructed spectral index maps. 
The lower most panel shows the difference between aver¬ 
age reconstructed index map and the input spectral index 

map, -/3j(p). Clearly, the difference map shows our 

method can accurately reconstruct spectral index map inside 
the pmlO mask. The reconstruction difference towards higher 
(lower) galactic coordinates becomes slightly larger due to 
weak polarized signal in these regions. 

Eor any given Monte-Carlo simulation x^Wos,l3d) values 
shows a clear minimum, Xmm^ some values of (Pos,/3d) = 
(^0s.min,/3rf,min)- Although a topographical plot of x^ifios^Pd) 
depends upon the particular random simulation under con¬ 
sideration, the features of such plots that arise on the av¬ 
erage due to foregrounds, CMB and detector noise, can be 
studied by making a plot of {x^iPos,/3d)) (where the av¬ 
erage is performed over all 1000 simulations for each pair 
WosjPd)) with respect to both indices. We have shown this 
at the top panel of Eigure The elongated contours of this 
plot shows the relatively higher sensitivity of the data with 
the variation of synchrotron spectral index compared to the 
same variation of dust spectral index. In the bottom panel 
of this hgure we show standard deviation map of x^ val¬ 
ues for all indices. Lower standard deviation values around 
{Xmin} indicates minimum is recovered satisfactorily by our 
method. Using results from Monte Carlo simulations we es¬ 
timate average value of (;Soi,min,/3rf,min) corresponding to xL« 
as (-2.84±0.01,2.00± 0.004). Small error on the spectral in- 
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Fig. 3.— Top panel: Variation of mean Pd) with respect to spectral 

indices obtained from 1000 Monte Carlo simulations. Bottom panel: Stan¬ 
dard deviation map of x^(/5o.n Pd) for any one of the 1000 simulations. 


dex values corresponding to Xm/n shows that our method can 
reliably recover global spectral indices. 

Using the Monte Carlo simulations we hnd that values of 
Xmm is distributed over a narrow range of 0.82 to 1.05. A 
plot of histogram of Xm,„ is shown in top panel of Eigure 
We estimate, (xLn) =0.93, averaging over 1000 Monte-Carlo 
simulations, with sample standard deviation 0.04. The proba¬ 
bility density function of this figure, follows a very symmetric 
distribution for xLn- This can be expected due to very large 
number of degrees of freedom (~ 4260) available in the analy¬ 
sis. In the bottom panel of this figure we have shown probabil¬ 
ity distribution of 20° smoothed Aside = 8 input spectral index 
map, from the region that survives after application of pmlO 
mask. A bimodal structure of the distribution is readily iden¬ 
tifiable corresponding to spectral indices ^ -2.85 and -3.16 
respectively. Comparing these values with the mean recon¬ 
structed spectral index -2.84 ±0.01, it is interesting to note 
that our method picks out the peak corresponding to a flatter 
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Ps{p) 


Fig. 4.— Top panel: Probability density of minimum values obtained 
from 1000 simulations. Red vertical line shows the location of sample mean 
corresponding to X™„ = 0.93. Light golden band around the mean value 
shows regions within Icr deviation. Bottom panel: Probability density func¬ 
tion of pixel indices obtained from smoothed input spectral index map after 
application of pmlO mask. 

index, from available, a pair of doubly degenerate ones. 

6.2. Weight 

One very important component of this work are the 
weights. We show weights for reconstruction of Qosip) and 
Qosip)APs{p) templates for the entire range of spectral in¬ 
dices, respectively in the top and bottom panel of Figure]^ At 
the low frequency side, where synchrotron emission is ^m- 
inant weights tend to be positive for all indices. We see this 
for K, A and Ka bands. For all other bands except F band, 
weights are negative for all indices. For F band weights be¬ 
come positive for all indices to cancel out signals from CMB, 
thermal dust, Qosip)^l3s{p) and at the same time to preserve 
Qos- For reconstruction of Qosip)APs{p) template weights 
must completely project out its shape vector from the data, 
at the same time, the weight vector need to be orthogonal 
to shape vectors of synchrotron and other components. Be¬ 
cause of close analytical similarity in the shape vectors of 
Qosip) and Qos(p)A/3s{p) the orthogonality condition leads to 
somewhat larger magnitude of weights for Qoi;{p)Ap^(p) than 
Qosip) template. Relatively larger magnitude of weights can 
easily be seen from the second panel than the hrst one. For 
the second panel and for K band, weights take largest neg¬ 
ative value, since, from the expression of shape vector for 


Qosip)Al3,{p) component, (a(iy)(iy/iyo,)~^0'\n{i^/iyos)), shape 
vector has a zero component at K band. For low frequency A, 
Ka, Q, B bands all weights are positive indicating its larger 
shape vector at lower frequency. For V, C, D and E bands 
weights become negative, since shape vector decreases in this 
frequency range. Finally for F band weights again become 
positive to project out QQs{p)Aj3s{p) component completely 
and to cancel out rest of the components. 

6.3. CMB and Foreground 

As associated results of our method, we reconstruct ML 
templates for all components using the best ht spectral index 
valu es I5 qs. mim Pd.min Corresponding to x»;m obtained in Sec- 
tion |6.1| We show the results in Figure]^ Unlike the case 
for reconstructed synchrotron spectral map, we do not apply 
any mask on the recovered templates, since CMB and fore¬ 
ground signals are expected to be stronger than the spectral 
index variations. Top, middle and bottom panels of this hgure 
show results for CMB, thermal dust and synchrotron template 
reconstruction respectively, obtained from a randomly cho¬ 
sen Monte Carlo simulation with seed = 100. All the color 
maps of this hgure are in fiK unit (thermodynamic for CMB, 
antenna for foreground templates). Following usual row, col¬ 
umn convention of a matrix, we represent any image of this 
hgure by the pair (i,J). (1,1) image represents reconstructed 
CMB map against the input CMB map for the same simula¬ 
tion, shown as (1,2) image. As one can easily visually iden¬ 
tify the large scale features of these two maps are very nearly 
identical. (1,3) image shows the difference between recon¬ 
structed and input CMB map and hence is dominated by resid¬ 
ual detector noise. The noisy nature of the difference map 
is also indicated by the symmetric color table plotted below 
(1,3) image. The (1,4) image represents the standard devia¬ 
tion map estimated from set of all (1,3) maps resulting from 
1000 Monte Carlo simulations of component separation algo¬ 
rithm. As expected, the standard deviation map clearly indi¬ 
cates a scan modulated noise pattern as applicable for WMAP 
observations. There exist some foreground induced errors on 
the galactic plane. Such errors can be reduced by extending 
our method to a second order perturbative analysis. 

The (2,1) image represents reconstructed thermal dust tem¬ 
plate at 217 GHz for seed = 100. The reconstructed template 
matches very well with the input thermal dust template at the 
same frequency, viz., (2,2) image. The difference between re¬ 
constructed and input template ((2,3) image) is dominated by 
noise. The standard deviation map for reconstructed thermal 
dust template is shown as (2,4) image. Clearly, WMAP’s scan 
modulated noise pattern is visible in this hgure. There also ex¬ 
ist some residual foregrounds error in this map, which again 
can be eliminated by performing a higher order perturbative 
analysis. 

The (3,1) image shows reconstructed synchrotron template 
at 23 GHz for the same random seed. The input synchrotron 
template at 23 GHz is shown as (3,2) image. Both these maps 
match with each each other very well. The difference of the 
hrst and second maps is shown as (3,3) image. The standard 
deviation map, image (3,4), shows the scan modulated de¬ 
tector noise pattern, without any visible signature of residual 
foreground error. 

7. DISCUSSION & CONCLUSION 

Reconstruction of CMB and foreground components jointly 
from microwave observations is a primary task for cosmolo- 
gists. This can be achieved relatively easily if i) observed 
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Fig. 5.— Top panel: Weights for reconstruction of Qos(p)- Bottom panel: Weights for reconstruction of Qo.i(p)^A(p). For any given panel, plots in consecutive 
rows from left to right indicate respectively weights for K, A, Ka, Q, B, C, V, D, E and F frequency bands. 


sky signal consists of only those components that have rigid 
frequency scaling all over the sky ii) detector noise is negli¬ 
gible and iii) total number of these components is less than 
or equal to total number of observed frequency bands. How¬ 
ever, the task is complicated when spectral index of any one 
of the components varies with sky positions. The stronger 
the variation effectively one gets more components with rigid 
frequency scaling, which, if not taken into account in the anal¬ 
ysis, may cause the results to be biased, or at the least difficult 
to interpret. The scenario becomes even worse in the presence 
of non-negligible detector noise. We have shown that in the 
presence of foregrounds following power law model with spa¬ 
tially varying index and small amount of detector noise one 
can recover the spectral index variation map reliably from the 
WMAP and Planck Stokes Q observations without any need 
to model foregrounds and or CMB in terms of any model tem¬ 
plates. The only assumption we have made in this work is a 


power law model for the foregrounds. However, this is not a 
necessary condition for our method to work, since one can di¬ 
rectly fit for the shape vectors for all components by expand¬ 
ing the parameter space. Then with the help of perturbative 
expansion of the spectral index variation one can reconstruct 
the index variation map. Since the spectral index variation 
for synchrotron component is small compared to other emis¬ 
sions (specifically, compared to synchrotron emission ampli¬ 
tude itself), and can be masked by the presence of compara¬ 
ble level of detector noise contamination, we chose to work 
with maps containing low level of detector noise, by smooth¬ 
ing the frequency maps by a 20° beam window function. We 
use full pixel-pixel detector noise covariance matrix to accu¬ 
rately model the smoothed detector noise pattern, and thus 
for better reconstruction of spectral index map. The condi¬ 
tion for large beam window can be relaxed easily for exper¬ 
iments, or for sky regions with negligible detector noise. A 
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Fig. 6.— Results for ML template reconstruction for CMB (top panel), thermal dust (middle panel) and synchrotron (bottom panel). Units for all colormaps 
are in temperature unit (thermodynamic for CMB, antenna for foreground components). 


very interesting future project will be to apply the method on 
the Planck polarized observations combined with WMAP, and 
with polarization sensitive, low noise, future CMB observa¬ 
tions like CMBPol. For very low noise experiments in future, 
it would be very useful to extend our method up-to second or¬ 
der perturbative analysis for even more accurate understand¬ 
ing of cosmological and astrophysical information contained 
in data. 

We summarize main results and conclusions of this work as 
follows: 


are ML estimates. 

4. The perturbative method of this work possesses an in¬ 
teresting property that it can be tuned to the data with an 
accuracy, as allowed by the data. This could be particu¬ 
larly very useful for CMB and foreground component 
separation problem if the data is insufficient to con¬ 
strain an otherwise completely accurate model, due to 
underlying degeneracies in the foreground components. 


1. We present and implement a first order perturbative 
method to estimate variation of synchrotron spectral in¬ 
dex over the microwave sky using simulated Stokes Q 
observations of WMAP K, Ka, Q, V, Planck LFI, and 
three lowest frequency Planck HFI frequency bands. 
We validate the method by performing 1000 Monte 
Carlo simulations over regions of the sky dominated by 
synchrotron Stokes Q signal. Our results suggest that, 
following the first order analysis synchrotron spectral 
index variation map can be reconstructed reliably by 
our method with a narrow error margin. 

2. As accompanying results of this work, we reconstruct 
CMB, thermal dust and synchrotron template compo¬ 
nents. An error estimation on the reconstructed tem¬ 
plates suggest that the template components can be re¬ 
constructed reliably by our method. 

3. Our reconstructed templates and best ht spectral indices 


The work presented in this paper will allow an improved es¬ 
timation of cosmological signal due to its detailed foreground 
modelling in terms of variation of spectral indices. It will be 
useful to extract astrophysical information about galactic syn¬ 
chrotron emission, in particular, information about the galac¬ 
tic cosmic ray electrons and the magnetic held. Finally, the 
method can be used to extract jointly all foregrounds compo¬ 
nents along with CMB signal from CMB observations. 
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crowave Background Data Analysis (LAMBDA), part of 
the High Energy Astrophysics Science Archive Center 
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trophysics Science Division at the NASA Goddard Space 
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